Fast track to GRASS with R: the rgrass package

Author

Veronica Andreo

The rgrass package allow us to interact with GRASS GIS tools (and data) serving as an interface between GRASS GIS and R. The rgrass package is developed and maintained by Roger Bivand and can be found at: https://github.com/rsbivand/rgrass/. In this fast track tutorial, we will learn how to use GRASS GIS from R.

Note

Setup note:

To run this tutorial locally you should have GRASS GIS 8.4+, R and, optionally, RStudio installed. You will also need to install rgrass, terra and mapview R packages and download the North Carolina sample dataset.

rgrass main funtions

The main functions within rgrass are the following:

initGRASS(): starts a GRASS GIS session from R.
execGRASS(): executes GRASS GIS commands from R.
gmeta(): prints GRASS GIS session metadata like database, project, mapset, computational region settings and CRS.
read_VECT() and read_RAST(): read vector and raster maps from a GRASS project into R.
write_VECT() and write_RAST(): write vector and raster objects from R into a GRASS GIS project.

Note

For further details on rgrass functionality, usage examples and data format coercion, see: https://rsbivand.github.io/rgrass/.

Basic Usage: Choose your own adventure

If you are a regular R user that needs to use GRASS GIS functionality because, well, you know it rocks, rgrass has your back. For example, maybe you struggle with large raster datasets in R or you need some specific tool, like watershed delineation for a large high resolution DEM. We will show here the way to use GRASS GIS tools within your R workflows.

On the other hand, if you already use GRASS as your geospatial data processing engine, you most likely have your spatial data within GRASS projects. You might need however to do some statistical analysis, some modelling and prediction or create publication ready visualizations in R. In such cases, you can start a GRASS GIS session in your project from R or RStudio.

Let’s see the general basic steps and then dive into the details:

Open R (or RStudio)
Load rgrass library with library(rgrass)
Start a GRASS GIS session with initGRASS()
Use GRASS GIS modules through execGRASS()
Use read_VECT(), read_RAST(), write_VECT() and write_RAST() to read data from and write data into GRASS database.

Note

GRASS raster and vector maps are translated into terra’s packege SpatRaster and SpatVector objects, respectively. These objects can then, within R, be easily coerced to other types of spatial objects such as simple features (sf), stars, etc.

See terra vignettes with further explanations and examples: https://rspatial.github.io/terra/.

A. Use GRASS GIS tools within your R spatial workflows

In the case you need to include some of the cool GRASS GIS tools within your R workflow, the initGRASS() function allows you to create temporary GRASS projects to use GRASS modules on R objects. This is equivalent to what QGIS does when you use GRASS tools via the Processing Toolbox.

So, the workflow goes like this: we will use initGRASS() to create a temporary GRASS project based on the extent, resolution and CRS of a raster or vector R object, likely the one we want to process or one that has the extent of our study area. Hence, we need to pass a reference spatial grid via the SG parameter. Then, we’ll write our R objects into the temporary GRASS project, run the desired processes, export the outputs back to the R environment and done.

Let’s start with getting some spatial data, eg. a raster file, into R.

library(terra)

terra 1.7.71

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
plot(r)

Now, we’ll load the rgrass library and start GRASS GIS with SG = r, so a GRASS project is internally created with r’s object CRS (BTW, you can check that with crs(r)), extent and resolution. These latter define the GRASS GIS computational region that will affect all raster processing, i.e., all new raster maps generated within GRASS GIS will have the same extent and resolution of the map provided. If you wish to change the computational region later on, you can use the g.region GRASS tool with execGRASS("g.region --h").

library(rgrass)

GRASS GIS interface loaded with GRASS version: (GRASS not running)

Note that we also pass tempdir() as the folder where our GRASS temporary project will be created. tempdir() will default to whichever is the temporary folder in your operative system.

# path to GRASS binaries
grassbin <- system("grass --config path", intern = TRUE)

# start GRASS GIS from R
initGRASS(gisBase = grassbin, 
          home = tempdir(),
          SG = r, 
          override = TRUE)

gisdbase    /tmp/Rtmp7Wq6Wd 
location    file5cda8663a8f01 
mapset      file5cda8325d3f4f 
rows        90 
columns     95 
north       50.19167 
south       49.44167 
west        5.741667 
east        6.533333 
nsres       0.008333333 
ewres       0.008333326 
projection:
 GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Now, we can write our SpatRaster into the GRASS GIS temporary project.

write_RAST(r, "terra_elev")

Importing raster map <terra_elev>...
   0%   3%   6%  10%  13%  16%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%

SpatRaster read into GRASS using r.in.gdal from file

Alternatively, we can use GRASS GIS importing tools to import common raster and vector formats. These will provide re-projection on the fly, if needed.

execGRASS("r.import", input=f, output="test")

Importing raster map <test>...
   0%   3%   6%  10%  13%  16%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%

Let’s check both are indeed within the project and run the GRASS module r.slope.aspect on it.

execGRASS("g.list", type = "raster")

terra_elev
test

execGRASS("r.slope.aspect", 
          elevation = "terra_elev", 
          slope = "slope",
          aspect = "aspect")

   0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  34%  37%  40%  43%  46%  49%  52%  56%  59%  62%  65%  68%  71%  74%  78%  81%  84%  87%  90%  93%  96% 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete

execGRASS("g.list", type = "raster")

aspect
slope
terra_elev
test

Let’s get slope and aspect maps into R

grass_maps <- read_RAST(c("aspect", "slope"))

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda866f11586.grd> created.
Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda858f9b728.grd> created.

grass_maps

class       : SpatRaster 
dimensions  : 90, 95, 2  (nrow, ncol, nlyr)
resolution  : 0.008333326, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : file5cda866f11586.grd  
              file5cda858f9b728.grd  
names       :       aspect,      slope 
min values  :   0.08974174, 0.01416342 
max values  : 360.00000000, 7.22943783

Now that the output maps are back into our R environment, we can plot them, do further analysis or write them into other raster formats, in which case we use terra::writeRaster() function.

plot(grass_maps)

writeRaster(grass_maps, "grass_maps.tif", overwrite=TRUE)

Alternatively, we can use GRASS GIS exporting tools like r.out.gdal and v.out.ogr, to directly save our outputs into common raster or vector files respectively.

execGRASS("r.out.gdal", input="slope", output="slope.tif", format="GTiff", flags="overwrite")

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
Using GDAL data type <Float32>
Input raster map contains cells with NULL-value (no-data). The value nan
will be used to represent no-data values in the input map. You can specify
a nodata value with the nodata option.
Exporting raster data to GTiff format...
ERROR 6: slope.tif, band 1: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
r.out.gdal complete. File <slope.tif> created.

B. Use R tools within GRASS GIS workflows

Let’s see an example for the case when we do all our geospatial data processing within GRASS and hence have all the spatial data organized within GRASS projects but we need to run some statistical analysis, modelling, prediction or visualization in R.

We start R or Rstudio and load the rgrass library. It will tell us that GRASS is not running, but we know that already… and that’s about to change in a moment.

library(rgrass)

So, we start GRASS from within R or RStudio using the initGRASS() function. Since we want to start GRASS in a specific project and mapset, we need to specify them. Optionally, we specify which GRASS binary to use. This might be useful in the case we have many GRASS versions in our system. If not provided, initGRASS() will attempt to find it in default locations depending on your operative system.

We can now list and read our GRASS raster and vector maps into R and do our statistical analysis, modelling and/or visualizations using other R packages. Obviously, we can also use rgrass as a mere interface to analyze data within GRASS, i.e., no reading from or writing to GRASS GIS needed, we can just use GRASS from R. Here, we’ll demonstrate the use of all the main rgrass functions mentioned above.

Let’s then list our GRASS raster and vector maps:

# list GRASS raster maps
execGRASS("g.list", parameters = list(type="raster"))

basins
developed
elevation
elevation_2
elevation_shade
geology
lakes
landuse
result_from_R
soils

# list GRASS vector maps
execGRASS("g.list", parameters = list(type="vector"))

boundary_region
boundary_state
census
elev_points
firestations
geology
geonames
hospitals
points_of_interest
railroads
roadsmajor
schools
streams
streets
zipcodes

The resulting map lists could be saved in an R object that we can subset later in case we want to import several but not all raster maps, for example. Let’s see how to do that.

# save map list in an object
rast_list <- execGRASS("g.list", parameters = list(type="raster"))

basins
developed
elevation
elevation_2
elevation_shade
geology
lakes
landuse
result_from_R
soils

rast_list

[1] 0
attr(,"resOut")
 [1] "basins"          "developed"       "elevation"       "elevation_2"    
 [5] "elevation_shade" "geology"         "lakes"           "landuse"        
 [9] "result_from_R"   "soils"          
attr(,"resErr")
character(0)

# retrieve only the maps list and overwrite rast_list
rast_list <- attributes(rast_list)$resOut

# import elevation and landuse
to_import <- c("elevation", "landuse") # optionally, by position: rast_list[c(3,7)]

maplist <- list()
for (i in to_import) {
  maplist[i] <- read_RAST(i)
}

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda8613bd54d.grd> created.

Warning in `[<-`(`*tmp*`, i, value = read_RAST(i)): implicit list embedding of
S4 objects is deprecated

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
Using GDAL data type <UInt32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda8a157673.grd> created.

Warning in `[<-`(`*tmp*`, i, value = read_RAST(i)): implicit list embedding of
S4 objects is deprecated

maplist

$elevation
class       : SpatRaster 
dimensions  : 1350, 1500, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(HARN) / North Carolina (EPSG:3358) 
source      : file5cda8613bd54d.grd 
name        : file5cda8613bd54d 
min value   :          55.57879 
max value   :         156.32986 

$landuse
class       : SpatRaster 
dimensions  : 1350, 1500, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(HARN) / North Carolina (EPSG:3358) 
source      : file5cda8a157673.grd 
name        : file5cda8a157673 
min value   :                0 
max value   :                7

Remember that raster objects will always be exported from GRASS GIS following the computational region settings. So, bear that in mind when reading into R which will hold them in memory. Vectors however will be exported in their full extent.

Let’s load the terra library to quickly display our recently imported raster maps:

library(terra)
plot(maplist$elevation)

Optionally, we could stack our two SpatRaster objects together and plot them together:

rstack <- rast(maplist)
plot(rstack)

Let’s create a boxplot of elevation per land class.

boxplot(rstack$elevation, rstack$landuse, maxcell=50000)

Let’s import a vector map, too, and explore its attributes.

census <- read_VECT("census")

Exporting 2547 areas (may take some time)...
   5%  11%  17%  23%  29%  35%  41%  47%  53%  59%  65%  71%  77%  83%  89%  95% 100%
WARNING: 10 features without category were skipped. Features without
         category are written only when -c flag is given.
v.out.ogr complete. 2537 features (Polygon type) written to <census> (GPKG
format).

head(census)

  cat OBJECTID      AREA PERIMETER BLOCK_ BLOCK_ID        BLOCKNUM TOTAL_POP
1  86    82802  9981.633   444.547  82803    82802 371830535104051        13
2  88    82805   304.076    91.371  82806    82805 371830535104046         0
3  91    82811  2400.879   192.304  82812    82811 371830535104034        20
4 140    82955  1239.619   147.011  82956    82955 371830535104039         0
5 173    83073 10358.715   429.738  83074    83073 371830535014007        80
6 228    83352 15459.671   458.961  83353    83352 371830535015007        18
  POP_1RACE WHITE_ONLY BLACK_ONLY AMIND_ONLY ASIAN_ONLY HWPAC_ONLY OTHER_ONLY
1        13         13          0          0          0          0          0
2         0          0          0          0          0          0          0
3        20          3          0          0         17          0          0
4         0          0          0          0          0          0          0
5        76         38          9          0          0          0         29
6        18         17          1          0          0          0          0
  POP_2RACES HISPANIC MALE FEMALE P_UNDER_5 P5_TO_9 P10_TO_14 P15_TO_19
1          0        0    6      7         1       2         1         1
2          0        0    0      0         0       0         0         0
3          0        0   10     10         1       1         1         1
4          0        0    0      0         0       0         0         0
5          4       54   37     43        20       4         5         7
6          0        3   10      8         1       0         1         0
  P20_TO_24 P25_TO_34 P35_TO_44 P45_TO_54 P55_TO_59 P60_TO_64 P65_TO_69
1         1         3         1         1         1         0         1
2         0         0         0         0         0         0         0
3         2         8         2         3         0         1         0
4         0         0         0         0         0         0         0
5        14        17         7         5         0         0         0
6         0         5         0         5         1         2         2
  P70_TO_74 P75_TO_84 P75_OLDER P85_OLDER MEDIAN_AGE P18_OLDER P21_OLDER
1         0         0         0         0       25.5         9         8
2         0         0         0         0        0.0         0         0
3         0         0         0         0       30.0        16        16
4         0         0         0         0        0.0         0         0
5         0         1         1         0       21.8        48        44
6         0         1         1         0       52.5        16        16
  P65_OLDER M75_OLDER F75_OLDER HOUSEHOLDS HH_SIZE SINGLE_HH FAMILY_HH
1         1         0         0          3    4.33         0         2
2         0         0         0          0    0.00         0         0
3         0         0         0          7    2.86         0         6
4         0         0         0          0    0.00         0         0
5         1         0         1         17    4.71         1        16
6         3         0         1          9    2.00         3         6
  MAR_COUPLE W_OWN_CHLD M_SNGL_PAR F_SNGL_PAR NONFAM_HH HH_W_CHILD HH_W_ELDER
1          2          2          0          0         1          3          1
2          0          0          0          0         0          0          0
3          6          3          0          0         1          3          0
4          0          0          0          0         0          0          0
5          5          3          3          3         0         10          1
6          6          2          0          0         0          2          2
  SNGL_ELDER POP_IN_HH FAM_HHLDER SPOUSE CHILD GROUP_QTRS GQ_INST GQ_NONINST
1          0        13          2      2     4          0       0          0
2          0         0          0      0     0          0       0          0
3          0        20          6      6     6          0       0          0
4          0         0          0      0     0          0       0          0
5          0        80         16      5    22          0       0          0
6          1        18          6      6     3          0       0          0
  FAM_SIZE HOUSING_U OCCUPIED VACANT OWNER_U RENTER_U VACANT_SEA HH_SIZE_OW
1     4.00         3        3      0       0        3          0       0.00
2     0.00         0        0      0       0        0          0       0.00
3     3.00         9        7      2       7        0          0       2.86
4     0.00         0        0      0       0        0          0       0.00
5     3.81        17       17      0       9        8          0       5.00
6     2.50         9        9      0       8        1          0       1.88
  HH_SIZE_RE
1          4
2          0
3          0
4          0
5          4
6          3

summary(census$TOTAL_POP)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.00   26.00   63.03   62.00 6173.00

plot(census, "P25_TO_34", type="interval", breaks=5, plg=list(x="topright"))

Let’s do some interactive visualization with mapview.

library(mapview)
mapview(rstack$elevation) + census

We highly recommend you to check the tmap package to make really appealing and publication ready maps.

To exemplify the use of write_* functions, let’s do a simple operation with the landuse raster map. We will apply a custom function that makes NULL all values less than 4.

result <- app(rstack$landuse, fun=function(x){ x[x < 4] <- NA; return(x)} )
plot(result)

write_RAST(result, "result_from_R", overwrite = TRUE)

WARNING: Raster map <result_from_R> already exists and will be overwritten
Importing raster map <result_from_R>...
   0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%

execGRASS("g.list", parameters = list(type="raster", pattern="result*"))

result_from_R

Finally, there’s yet another way in which you can use GRASS and R together, and it involves calling R from the GRASS terminal. In this way, rgrass will read all GRASS GIS session environmental variables, and you won’t need to use initGRASS. It goes more or less like this:

Open GRASS GIS
Type R or rstudio & in the GRASS GIS terminal
Load rgrass library with library(rgrass)
Use read_VECT(), read_RAST() to read data from GRASS into R
Do your analysis or plotting in R
Write data (back) to GRASS database with write_VECT() and write_RAST()
Quit R quit() and get back to GRASS GIS terminal.

Starting GRASS GIS...

          __________  ___   __________    _______________
         / ____/ __ \/   | / ___/ ___/   / ____/  _/ ___/
        / / __/ /_/ / /| | \__ \\_  \   / / __ / / \__ \
       / /_/ / _, _/ ___ |___/ /__/ /  / /_/ // / ___/ /
       \____/_/ |_/_/  |_/____/____/   \____/___//____/

Welcome to GRASS GIS 8.3.0
GRASS GIS homepage:                      https://grass.osgeo.org
This version running through:            Bash Shell (/bin/bash)
Help is available with the command:      g.manual -i
See the licence terms with:              g.version -c
See citation options with:               g.version -x
If required, restart the GUI with:       g.gui wxpython
When ready to quit enter:                exit

Launching <wxpython> GUI in the background, please wait...
[Raster MASK present]
GRASS nc_basic_spm_grass7/PERMANENT:~ > R

R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(rgrass)
GRASS GIS interface loaded with GRASS version: GRASS 8.3.0 (2023)
and location: nc_basic_spm_grass7
>

References

Bivand R (2024). rgrass: Interface Between ‘GRASS’ Geographical Information System and ‘R’. R package version 0.4-1, https://CRAN.R-project.org/package=rgrass.